213-Discussion

Discussion I

For today’s discussion, we will review the distinction between parametric and nonparametric regression models, as well as some key properties of simple linear regression. If you wish to review materials from POL 211 & 212, you can access them through the following links: POL 211 Discussion, POL 212 Discussion.

1. Parametric versus nonparametric Regression Models
2. Simple Linear Regression
3. Properties of Least Squares Estimator

1. Parametric versus nonparametric Regression Models

When it comes to regression analysis, choosing the right approach is crucial for accurate predictions and meaningful insights. Two common methods used are parametric, like linear regression, and semi/non-parametric, like smoothing spline regression or Kernal regression. Each has its own advantages and disadvantages, and the choice between them largely depends on the nature of the data and the underlying relationships.

a. What is parametric and semi/nonparametric regression?

Parametric Regression Linear regression is a well-known parametric method that assumes a linear functional form for the relationship between the predictors (\(X\)) and the target variable (\(Y\)). This approach has several benefits, such as ease of estimation with a small number of coefficients. In linear regression, these coefficients have straightforward interpretations, and statistical significance tests are readily applicable. However, parametric methods come with a significant limitation — they rely on the assumption that the specified functional form is a close approximation to the true relationship. If this assumption is far from reality, linear regression can perform poorly and yield unreliable results.

Nonparametric Regression On the other hand, non-parametric methods like K-Nearest Neighbors (KNN) regression do not make explicit assumptions about the functional form of the relationship between \(X\) and \(Y\). Instead, they provide a more flexible approach for regression. KNN regression identifies the K training observations closest to a prediction point and estimates the target variable by averaging their values. While this approach is more versatile and can handle complex relationships, it can suffer from high variance when K is small, leading to overfitting. Conversely, when K is large, KNN regression can underfit the data.

Example

Assume that we have a outcome variable \(Y\) and two explanatory variables, \(x_1\) and \(x_2\). In general, the regression model that describes the relationship can be written as:

\[Y = f_1(x_1) + f_2(x_2) + \epsilon\]

• Some parametric regression models:

  • Multiple linear regression model: \(Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon\)
  • Polynomial regression model of second order: \(Y = \beta_0 + \beta_{10}x_1 + \beta_{11}{x_1}^2 + \beta_{20}x_2 + \beta_{21}{x_2}^2 + \epsilon\)
  • Nonlinear regression model: \(Y = \beta_0 + \beta_1x_1 + \beta_2e^{\beta_3x_2} + \epsilon\)
  • Poisson regression with \(Y\) is count: \(log(\mu) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \epsilon\)

• If we do not know \(f_1\) and \(f_2\) functions, we need to use a Nonparametric regression model.

b. Nonparametric regression estimation methods: K-Nearest Neighbors (KNN) regression

K-Nearest Neighbors (KNN) regression is one of the simplest and best-known nonparametric methods. Given a value for \(K\) and a prediction point \(x_0\), KNN regression first identifies the \(K\) training observations that are closest to \(x_0\), represented by \(N_0\). It then estimates \(f(x_0)\) using the average of all the training responses in \(N_0\). In other words,

\[\hat{f}(x_0)=\frac{1}{K}\sum_{x_i \in N_0}y_i\]

When using KNN as a regressor with a continuous dependent variable, data points are scattered across the coordinate plane. When a new data point is introduced, the number of neighbors (\(K\)) is determined using any of the distance metrics.Usually, the Euclidean distance is used as the distance metric. Once the neighbors are identified, the predicted value of the new data point is calculated as the average of all the neighbors’ values combined.

The below figure illustrates two KNN fits on a data set with \(p = 2\) predictors. The fit with \(K = 1\) is shown in the left-hand panel, while the right-hand panel corresponds to \(K = 9\).We see that when \(K = 1\), the KNN fit perfectly interpolates the training observations, and consequently takes the form of a step function. When \(K = 9\), the KNN fit still is a step function, but averaging over nine observations results in much smaller regions of constant prediction, and consequently a smoother fit. In general, the optimal value for \(K\) will depend on the bias-variance tradeoff, which Chris introduced in POL 212. A small value for \(K\) provides the most flexible fit, which will have low bias but high variance. This variance is due to the fact that the prediction in a given region is entirely dependent on just one observation. In contrast, larger values of \(K\) provide a smoother and less variable fit; the prediction in a region is an average of several points, and so changing one observation has a smaller effect. However, the smoothing may cause bias by masking some of the structure in \(f(X)\).

The bias-variance tradeoff is a fundamental concept in machine learning that deals with the balance between the bias (error from overly simplistic assumptions) and variance (sensitivity to fluctuations in the training data) of a model. In the context of KNN (K-Nearest Neighbors), adjusting the value of K can help control this tradeoff.

• Low \(K\) (e.g., \(K=1\)): Low bias, high variance. The model closely follows the training data, leading to high variance and sensitivity to noise. It can capture complex patterns but may overfit the training data.

• High \(K\) (e.g., large \(K\)): High bias, low variance. The model averages over more data points, resulting in smoother predictions with lower variance but potentially higher bias. It may underfit the data by oversimplifying the underlying patterns.

Choosing an appropriate value for \(K\) involves finding a balance between bias and variance to achieve optimal model performance. This tradeoff influences the model’s ability to generalize to unseen data.

There are no pre-defined statistical methods to find the most favorable value of \(K\), but here is a typical way of how to choose an optimal \(K\):

1. Initialize a random \(K\) value and start computing.
2. Choosing a small value of \(K\) leads to unstable decision boundaries.
3. The substantial \(K\) value is better for classification as it leads to smoothening the decision boundaries.
4. Derive a plot between error rate and \(K\) denoting values in a defined range. Then choose the \(K\) value as having a minimum error rate.

c. When should we use parametric or nonparametric regression?

The key question is when to choose a parametric approach like linear regression over a non-parametric one such as KNN regression. The answer is straightforward: a parametric approach performs better when the chosen functional form is a close match to the true relationship, particularly in the presence of a linear relationship. If the specified functional form is far from the truth, and prediction accuracy is our goal, then the parametric method will perform poorly. For instance, if we assume a linear relationship between \(X\) and \(Y\) but the true relationship is far from linear, then the resulting model will provide a poor fit to the data, and any conclusions drawn from it will be suspect.

In contrast, non-parametric methods do not explicitly assume a parametric form for \(f(X)\), and thereby provide an alternative and more flexible approach for performing regression.

To illustrate this point, let’s consider a few scenarios and use R to simulate them:

1. Linear Relationship: When the true relationship between \(X\) and \(Y\) is linear, linear regression outperforms nonparametric regression. Linear regression provides an almost perfect fit in this situation, as it closely matches the underlying relationship.

   set.seed(1234) # for replication
   N <- 1000 # set up the sample size
   x1 <- rnorm(N, mean = 10, sd = 3)
   e <- rnorm(N, mean = 0, sd = 3) # set the errors to be normally distributed 
   y <- 7 + 3*x1 + e # set the true y
   df <- data.frame(id = 1:N, y = y, x1 = x1, e = e)

   In the above chunk, I specify the true \(y\) to be a function of \(x_1\) and \(e\). The true values of the regression coefficients is 3 for \(x_1\). The intercept is 7. Let’s first plot the true relationship between \(x_1\) and \(y\)!

   library(tidyverse)
   ggplot(data = df, aes(y = y, x = x1)) +
     geom_point(alpha = 0.6) +
     geom_smooth() +
     theme_bw()

   

   Let’s begin by splitting the generated data into a training set, comprising 80% of the data, and a test set, comprising the remaining 20%.

   set.seed(1234) # set the seed to make the partition reproducible
   train <- df %>% sample_frac(0.8) # should have 800 obs
   test <- anti_join(df, train, by = 'id') # should have 200 obs

   Now, let’s first apply the lm() function to our training set. We’ll then use the coefficients estimated from this model to predict the y-values in the test set, followed by computing the mean squared error (MSE).

   fit <- lm(y ~ x1, data = train)

   	OLS
   Variables	Estimates	S.E.	C.I. (95%)
   (Intercept)	6.50 ***	0.35	5.81 – 7.19
   x1	3.06 ***	0.03	2.99 – 3.13
   Observations	800
   R2 / R2 adjusted	0.910 / 0.910
   p<0.05   ** p<0.01   *** p<0.001

   test$y_hat <- predict(fit, newdata = data.frame(x1 = test$x1))
   
   # compute mean squared error (MSE)
   MSE_ols <- sum((test$y - test$y_hat)^2)/nrow(test)
   
   MSE_ols

   ## [1] 8.770813

   The mean squared error (MSE) obtained by using OLS regression to estimate the relationship between \(y\) and \(x_1\) in the training set and then using the estimated coefficients to predict \(y\) in the test set is 8.770813. Now, let’s proceed to use KNN regression with various values of \(K\) to observe how the MSE varies.

   We will use the knn.reg() function from the FNN package for KNN regression. This function requires four arguments to be specified:

   • train: the predictors of the training data
   • test: the predictor values, \(x\), at which we would like to make predictions
   • y: the response for the training data
   • k: the number of neighbors to consider

   # install.packages("FNN")
   library(FNN)
   train_no_y <- train %>% select(x1) 
   test_no_y <- test %>% select(x1) 
   train_y <- train %>% select(y) 
   
   yhat_k1 <- knn.reg(train = train_no_y, test = test_no_y, y = train_y, k = 1)
   yhat_k1

   ## Prediction:
   ##   [1] 29.381821 28.304426 65.609505 23.837209 34.330518 33.436274 21.824014
   ##   [8] 30.847990 32.289959 42.276634 48.391432 53.703199 31.595459 33.423440
   ##  [15] 32.935792 41.304789 26.486024 39.445252 52.119807 38.608137 44.000824
   ##  [22] 38.608137 37.159409 38.540173 33.079801 39.907787 44.810978 37.870818
   ##  [29] 34.348167 42.174368 40.401263 34.996959 27.478094 48.391432 33.978417
   ##  [36] 26.367731 56.238177 31.240796 37.950729 55.185947 34.850824 43.470660
   ##  [43] 44.483329 21.700094 49.108239 44.558484 41.288458 40.623165 31.965684
   ##  [50] 28.174148 35.309716 44.393660 40.512042 48.391432 26.800229 31.877765
   ##  [57] 45.551430 28.354259 39.391663 43.100323 33.168206 25.652004 37.857259
   ##  [64] 42.280182 31.905242 30.070676 45.434786 32.597936 55.713238 30.245094
   ##  [71] 35.747199 24.318503 35.309716 42.142535 44.687740 34.463743 41.607736
   ##  [78] 11.887226 36.074342 55.185947 33.106591 64.709783 20.566998 54.842595
   ##  [85] 45.434786 51.245243 39.073328 28.874519 41.177921 35.357025 30.071689
   ##  [92] 36.168858 35.095983 28.392045 31.241205 32.798233 29.709591 37.191013
   ##  [99] 38.608137 26.841422 44.442204 24.630053 26.367731 19.105273 33.290191
   ## [106] 40.401263 56.238177 21.708475 33.771862 38.654251 42.868435 49.161958
   ## [113] 33.168206 26.492572 27.573812 41.597145 30.479421 45.411737 40.623165
   ## [120] 37.208831 35.170920  8.539242 48.085100 44.483329 20.591604 21.824014
   ## [127] 45.068575 40.911258 35.964535 44.731334 38.432659 44.912912 38.520474
   ## [134] 33.390623 24.145923 22.517012 37.950729 41.465057 36.782189 31.241205
   ## [141] 33.423440 29.709591 28.827443 21.824014 39.739236 49.108239 41.337496
   ## [148] 51.204406 31.801673 31.406837 30.414360 32.690670 39.641522 28.351906
   ## [155] 44.206919 32.373809 48.402430 31.595459 67.724169 27.404886 36.782189
   ## [162] 42.020108 54.160505 36.074342 48.649611 34.330518 44.253326 42.549677
   ## [169] 36.377909 46.279074 27.157538 25.692918 30.150333 29.936407 50.422834
   ## [176] 27.404886 37.870818 33.978417 29.820616 32.646670 28.464324 26.137602
   ## [183] 29.774837 38.787417 29.728468 39.795625 54.842595 48.649611 38.169491
   ## [190] 38.146542 33.233128 36.284027 41.820483 35.838288 38.441584 32.682236
   ## [197] 28.521211 41.054304 52.592726 44.479609

   As observed from the above results, when setting \(K = 1\) in knn.reg(), it provides the predicted \(y\) values for the test set. Now, let’s calculate the MSE for this prediction.

   MSE_k1 <- sum((test$y - yhat_k1$pred)^2)/nrow(test)
   MSE_k1

   ## [1] 20.25002

   As you can see, the MSE when employing KNN regression with \(K = 1\) is significantly higher than the MSE observed with OLS regression. Let’s create a for loop and repeat this process for various values of \(K\) to determine whether the MSE improves or worsens compared to the MSE when using OLS.

   MSE_KNN <- c()
   
   for (i in seq(1, 100, by = 2)){
     yhat_k <- knn.reg(train = train_no_y, test = test_no_y, y = train_y, k = i)
     MSE_k <- sum((test$y - yhat_k$pred)^2)/nrow(test)
     MSE_KNN <- c(MSE_KNN, MSE_k)
   }
   
   MSE_KNN

   ##  [1] 20.250016 11.554227 10.170821 10.020766 10.002982  9.655185  9.468217
   ##  [8]  9.151752  9.058586  8.825635  8.704344  8.733093  8.672769  8.708741
   ## [15]  8.678995  8.627246  8.531880  8.596169  8.653730  8.691982  8.783636
   ## [22]  8.831042  8.858023  8.874769  8.833497  8.853275  8.927840  8.982527
   ## [29]  8.938121  8.948997  9.004996  9.012363  8.970862  8.969159  9.031482
   ## [36]  9.058308  9.110250  9.158928  9.255961  9.305864  9.352528  9.415657
   ## [43]  9.468395  9.539323  9.576483  9.650497  9.693073  9.724640  9.752314
   ## [50]  9.824734

   ggplot(data = data.frame(K = seq(1, 100, by = 2), MSE = MSE_KNN), 
      aes(y = MSE, x = K)) +
     geom_line(alpha = 0.4, color = "darkgreen") +
     geom_point(alpha = 0.3, color = "darkgreen") +
     geom_hline(yintercept = MSE_ols, color = "darkgray", 
            alpha = 0.8, linetype = "dashed") +
     annotate("text", x = 90, y = 18, label = "KNN", color = "darkgreen", alpha = 0.5) +
     annotate("text", x = 90, y = 17, label = "OLS", color = "darkgray") +
     theme_bw()

   

   After simulation, it’s evident that when the true relationship between the variables of interest is linear, KNN regression doesn’t offer superior predictions compared to OLS regression.

2. Slight Non-Linearity: In cases of slight non-linearity, where the true relationship deviates slightly from linearity, KNN regression can perform nearly as well as linear regression. It still provides reasonable results without a substantial reduction in prediction accuracy.

3. Strong Non-Linearity: However, in situations with a strong non-linear relationship, KNN regression outperforms linear regression. This is because KNN can adapt to complex relationships, providing more accurate predictions.

   set.seed(1234) # for replication
   N <- 1000 # set up the sample size
   x1 <- rnorm(N, mean = 1, sd = 10)
   e <- rnorm(N, mean = 3, sd = 10) # set the errors to be normally distributed 
   y <- 3 + 4*(x1+3)^2*(x1-1)^3 + e # set the true y
   df <- data.frame(id = 1:N, y = y, x1 = x1, e = e)

   For example, let’s now generate a relationship between \(x_1\) and \(y\) that is strongly non-linear. I specify the true \(y\) to be a function of \(x_1\) and \(e\). Overall, the formula, \(y = 3+4(x_1+3)^2(x_1-1)^3+e\), describes a relationship where \(y\) is influenced by \(x_1\) in a nonlinear way, with two distinct “bumps” or curves, along with some random variability represented by \(e\).

   library(tidyverse)
   ggplot(data = df, aes(y = y, x = x1)) +
     geom_point(alpha = 0.6) +
     geom_smooth() +
     theme_bw()

   

   Following our previous approach, let’s now compare the MSE obtained from OLS regression, used to estimate the relationship between \(y\) and \(x_1\), with the MSE derived from KNN regression, employing different values of \(K\).

   • Step 1: Split the Sample in to Training and Test Set

     set.seed(1234) # set the seed to make the partition reproducible
     train <- df %>% sample_frac(0.8) # should have 800 obs
     test <- anti_join(df, train, by = 'id') # should have 200 obs
     fit2 <- lm(y ~ x1, data = train)

   • Step 2: Run OLS Regression and Compute the MSE

     fit2 <- lm(y ~ x1, data = train)

     	OLS
     Variables	Estimates	S.E.	C.I. (95%)
     (Intercept)	95970.43	389692.62	-668973.27 – 860914.13
     x1	734839.02 ***	38319.47	659620.15 – 810057.88
     Observations	800
     R2 / R2 adjusted	0.315 / 0.315
     p<0.05   ** p<0.01   *** p<0.001

     test$y_hat <- predict(fit, newdata = data.frame(x1 = test$x1))
     MSE_ols <- sum((test$y - test$y_hat)^2)/nrow(test)
     
     MSE_ols

     ## [1] 4.742095e+13

   • Step 3: Run KNN and Compute the MSE

     train_no_y <- train %>% select(x1) 
     test_no_y <- test %>% select(x1) 
     train_y <- train %>% select(y) 

     MSE_KNN <- c()
     
     for (i in seq(1, 100, by = 2)){
       yhat_k <- knn.reg(train = train_no_y, test = test_no_y, y = train_y, k = i)
       MSE_k <- sum((test$y - yhat_k$pred)^2)/nrow(test)
       MSE_KNN <- c(MSE_KNN, MSE_k)
     }
     
     MSE_KNN

     ##  [1] 1.006088e+11 1.791096e+11 1.840990e+11 5.863832e+11 9.814868e+11
     ##  [6] 4.548076e+11 1.081960e+12 1.129925e+12 1.487238e+12 2.374122e+12
     ## [11] 2.763030e+12 3.688595e+12 4.638701e+12 5.402272e+12 6.302316e+12
     ## [16] 6.098534e+12 6.924194e+12 7.720989e+12 8.508775e+12 9.279200e+12
     ## [21] 9.522820e+12 1.021639e+13 1.021119e+13 1.086252e+13 1.149660e+13
     ## [26] 1.209597e+13 1.238768e+13 1.296902e+13 1.354249e+13 1.352895e+13
     ## [31] 1.408380e+13 1.458944e+13 1.475359e+13 1.524066e+13 1.571534e+13
     ## [36] 1.618749e+13 1.664102e+13 1.702627e+13 1.745704e+13 1.775741e+13
     ## [41] 1.815320e+13 1.818130e+13 1.853354e+13 1.892895e+13 1.930142e+13
     ## [46] 1.966064e+13 2.001759e+13 2.037685e+13 2.072546e+13 2.106549e+13

   • Step 4: Compare the MSE obtained from OLS and with that from KNN

     ggplot(data = data.frame(K = seq(1, 100, by = 2), MSE = MSE_KNN), 
            aes(y = MSE, x = K)) +
       geom_line(alpha = 0.4, color = "darkgreen") +
       geom_point(alpha = 0.3, color = "darkgreen") +
       geom_hline(yintercept = MSE_ols, color = "darkgray", 
          alpha = 0.8, linetype = "dashed") +
       annotate("text", x = 90, y = 4e+13, label = "KNN", color = "darkgreen", alpha = 0.5) +
       annotate("text", x = 90, y = 3.7e+13, label = "OLS", color = "darkgray") +
       theme_bw()

     

     As observed, regardless of the choice of \(K\) in our KNN regression, the MSE obtained from KNN is consistently much smaller than that from OLS.

4. Curse of Dimensionality: When dealing with high-dimensional data (i.e., when you have a lot of predictors), KNN regression may suffer from the “curse of dimensionality.” In such cases, the performance of KNN deteriorates significantly as the dimensionality of the data increases. In the case of KNN regression, as the number of dimensions increases, the distance between data points becomes less meaningful. This is because in high-dimensional spaces, the concept of distance becomes less discriminating, as most data points are located far apart from each other. Consequently, KNN may struggle to find relevant neighbors, leading to poor predictive performance and increased computational complexity. On the other had, linear regression, with its fewer parameters, is less affected by this issue.

   

   Even in problems in which the dimension is small, we might prefer linear regression to KNN from an interpretability standpoint. If the test MSE of KNN is only slightly lower than that of linear regression, we might be willing to forego a little bit of prediction accuracy for the sake of a simple model that can be described in terms of just a few coefficients, and for which p-values are available.

2. Simple Linear Regression

Below is a population regression line equation, assuming that we are the god and know the data-generating process of the relationship between \(X\) (independent/explanatory variable) and \(Y\) (dependent/response/outcome variable).

However, since we usually don’t have the population data but only have access to a sample of it, we can only use the sample dataset to estimate the population regression line. In the below figure, \(b_0\) is the estimate of the regression intercept \(\beta_0\), while \(b_1\) is the estimate of the regression coefficient/slope \(\beta_1\).

The Ordinary Least Squares (OLS) approach aims to fit a line that minimizes the sum of the squared residuals, which means that the goal is to find a pair of \(b_0\) and \(b_1\) that can minimize the difference between our observed \(Y_i\) and estimated \(\hat{Y_i}\).

\[S(A, B) = min \Sigma E_i = min \Sigma (Y_i-\hat{Y_i})^2=min \Sigma (Y_i-(A+BX_i))^2\]

\[S(A, B) = \Sigma (Y_i-A-BX_i)^2\]

To find a pair of \(b_0\) and \(b_1\) that can minimize the difference between our observed \(Y_i\) and estimated \(\hat{Y_i}\) is an optimization problem. All we have to do is to take the partial derivatives of the above equation with respect to \(b_0\) and \(b_1\), respectively, set them equal to zero, and then solve it.

\[\frac{\partial S(A, B)}{\partial A}=\Sigma (-1)(2)(Y_i-A-BX_i)=0\]

\[\frac{\partial S(A, B)}{\partial B}=\Sigma (-X_i)(2)(Y_i-A-BX_i)=0\]

By solving the above two equations, you will get:

\[A = \bar{Y}-B\bar{X}\]

\[B = \frac{\Sigma (X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma (X_i-\bar{X})^2}\]

a. Gauss-Markov Theorem

Gauss Markov theorem states that “if the errors are independently distributed with zero expectation and constant variance, then the least-squares estimator b is the most efficient linear unbiased estimator of \(\beta\). That is, of all unbiased estimators that are linear functions of the observations, the least-squares estimator has the smallest sampling variance and, hence, the smallest mean-squared error. For this reason, the least-squares estimator is sometimes termed BLUE, an acronym for best linear unbiased estimator” (Fox 2016: 212).

In other words, under certain conditions, the ordinary least squares (OLS) estimates are thought to be the best linear unbiased estimator (BLUE). Here, I review a set of assumptions (known as Gauss-Markov assumptions) that leads to a BLUE estimator.

1. Linearity: the expected (mean) value of the disturbance term is 0 (i.e., \(E(u_i)=0\)).

   Fox (2016) explains that a violation of this assumption implies “that the model fails to capture the systematic pattern of relationship between the response and explanatory variables” (307). This assumption could be broken in two ways: one, the mean of error could be a constant across all values of \(X_i\), two, it could vary. In the former case, the intercept term will be consistently off by the error amount. A uniform measurement error across all observations would be a cause of such an error (Berry 1993: 42). A common cause of the latter case will be omitted variable bias. Omitted variable bias deserves particular attention because it can cause the OLS estimate to be both biased and inconsistent. Another common cause is the wrong functional form (e.g., specifying a linear relation instead of a quadratic relation). Thus, when nonlinearity is detected, we should explore not only variable transformations and different model specifications (Fox, 2016: 317), but also the possibility of omitted variables. In short, before running OLS, reseachers should make sure that the inclusion of relevant variables in the model and the setting of the function are resonable and are based on the theoretical understanding.

2. Nonstochastic regressors: \(X\) values are independent of the error term (i.e., \(cov(X_i, u_i)=0\)).

   The violation of this assumption suggests that there is an endogeneity problem between the explanatory and outcome variables, which can arise from measurement error on \(X\), omitted confounders, or reciprocal causation between \(X\) and \(Y\). Again, if the main goal of the research is to provide an unbiased and consistent estimator to explain how a social phenomenon works in a real-life setting, to provide a convincing explanation requires us to make sure that the effect of our explanatory variable is actually exogenous.Fortunately, scholars have developed many methods to address this issue. For example, researchers could employ an instrumental variable or matching to improve the unbiasedness of the estimated effect.

3. Homoskedasticity: constant error variance across values of \(X_i\) (i.e., \(Var(\varepsilon)=\sigma^2\)).

   One instance where heteroskedasticity is often violated is cross-sectional studies (Berry, 1993: 73). For example, because more developed countries can have better data accuracy, measurement error could be correlated with the level of development. Thus, in a study that includes the level of development as one of the independent variables, the variance of error term may not be constant (Berry, 1993: 73). In such circumstances, Fox (2016) suggests that we can run weighted-least-squares (WLS) regression to overcome such probelm (304). Other solutions include transformation of the dependent variable and correction of coefficient standard errors (e.g. robust standard errors and clustered standard errors) for heteroskedasticity (Fox, 2016: 307).

4. Independence: the observations are sampled independently. no autocorrelation between disturbances (i.e., \(cov(u_i, u_j)=0\) for \(i \neq j\)).

   The independence assumption is important for time-series analysis, where the assumption of spherical errors is often violated. To overcome this, Beck & Katz (1995) propose panel-correlated standard errors to when analyzing time-series data.

5. Normality: errors are distributed normally (i.e.,\(\varepsilon \sim N(0, \sigma^2)\)).

Assumptions of 1 and 2 are related to the bias of the estimator, whereas assumptions of 3 and 4 are related to the efficiency of the estimator: when the assumption of 1 or 2 is violated, the OLS estimator is no longer unbiased; when the assumption of 3 or 4 is violated, the OLS estimator is no longer efficient. With these weak set of assumptions (1-4), according to the Gauss-Markov theorem, the OLS estimator is BLUE.

(In statistics, efficiency is a measure of the quality of an estimator, which can be characterized by having a smaller possible variance.)

With the fifth assumption, normality, the OLS estimator is the most efficient among all the unbiased estimators, not just linear estimators. Adding this fifth assumption makes the Gauss-Markov theorem a strong set. This last assumption is important when we aim to make a causal inference using the OLS estimator. With the normal distribution of the errors assumption, we can also infer that the OLS estimators also have normal sampling distributions. As a result, it allows us to apply statistical tests such as t-tests even for small sample size studies.

One important caveat is that when there is perfect multicollinearity (i.e. when two or more independent variables are perfectly correlated), we cannot even get an estimate using the OLS method. However, even when there is high multicollinearity, the OLS estimate is still BLUE (Berry 1993: 27). In such cases, standard errors will be very high, making the estimates fluctuate considerably from sample to sample (Berry 1993: 27).

The below table presents the consequences of the violation of Gauss-Markov assumptions and corresponded suggested solutions.

Assumption	Violation	Solution
Linearity	Biased/inconsistent estimates	Transformations, polynomials, different model
Nonstochastic regressors	Biased/inconsistent estimates	Instrumental variables
Homoskedasticity	Biased standard errors	Robust standard errors
Independence	Biased standard errors	Fixed effects
Perfect collinearity	cannot run OLS	Omit one collinear term

b. Properties of Least Squares Estimator

The sample least squares coefficients are unbiased estimators of the population regression coefficients.

Here are some reviews before demonstrating the unbiasedness of the least-squares estimators \(A\) and \(B\) for \(\alpha\) and \(\beta\):

• The assumed model: \(Y_i=\alpha+\beta X_i+\varepsilon_i\)
• Linearity assumption: \(E(\varepsilon_i)=0\)
• Normality assumption: \(\varepsilon_i\sim N(0, \sigma_\varepsilon^2)\)
• Independence assumption: \(Cov(\varepsilon_i, \varepsilon_j)=0, for i\neq j\)
• The expectation of linear combination: \(E(a+bY)=a+bE(Y)\)
• The variance of linear combination:
  • \(Var(a+bY)=b^2Var(Y)\)
  • \(Var(X+Y)=Var(X)+Var(Y)+2Cov(X, Y)\)
• For the below demonstration, the \(X_i\) are assumed to be fixed, not random.
  • \(B=\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}=\frac{\Sigma(X_i-\bar{X})Y_i-\Sigma(X_i-\bar{X})\bar{Y}}{\Sigma(X_i-\bar{X})^2}=\frac{\Sigma(X_i-\bar{X})Y_i-\bar{Y}\Sigma(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}=\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2}\)
  • \(\Sigma(X_i-\bar{X})^2=\Sigma(X_i-\bar{X})(X_i-\bar{X})=\Sigma X_i(X_i-\bar{X})-\bar{X}\Sigma(X_i-\bar{X})=\Sigma X_i(X_i-\bar{X})\)

1. Demonstrate the unbiasedness of the least-squares estimators \(B\) for \(\beta\) in simple regression (i.e., \(E(B) = \beta\)).

   \(E(B)=\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}\)

   \(=\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2}\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}E(\Sigma(X_i-\bar{X})Y_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})E(Y_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})E(\alpha + \beta X_i+\varepsilon_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})(\alpha +\beta X_i+E(\varepsilon_i))\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}\Sigma(X_i-\bar{X})(\alpha+\beta X_i)\)

   \(=\frac{1}{\Sigma(X_i-\bar{X})^2}(\Sigma(X_i-\bar{X})\alpha+\Sigma(X_i-\bar{X})\beta X_i)\)

   \(=\frac{\beta}{\Sigma(X_i-\bar{X})^2}(\Sigma(X_i-\bar{X})X_i)\)

   \(=\beta\)

2. Demonstrate the unbiasedness of the least-squares estimators \(A\) for \(\alpha\) in simple regression (i.e., \(E(A) = \alpha\)).

   \(E(A)=E(\bar{Y}-B\bar{X})\)

   \(=E(\frac{1}{n}\Sigma Y_i-\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}\bar{X})\)

   \(=E(\frac{1}{n}\Sigma Y_i-\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2}\bar{X})\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2})Y_i)\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2})(\alpha +\beta X_i +\varepsilon_i))\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}))E(\alpha +\beta X_i +\varepsilon_i)\)

   \(=E(\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}))(\alpha+\beta X_i+E(\varepsilon_i))\)

   \(=\Sigma(\frac{1}{n}-\frac{\bar{X}(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2})(\alpha+\beta X_i)\)

   \(=\Sigma \frac{\alpha}{n}+\Sigma \beta \frac{X_i}{n}-\alpha \frac{\bar{X}\Sigma(X_i-\bar{X})}{\Sigma(X_i-\bar{X})^2}-\beta \frac{\bar{X}\Sigma(X_i-\bar{X})X_i}{\Sigma(X_i-\bar{X})^2}\)

   \(=\alpha + \beta \bar{X}-\alpha*0-\beta \bar{X}\)

   \(=\alpha\)

3. Derive the sampling variances of \(A\) in a simple regression.

   \(Var(B)=Var(\frac{\Sigma(X_i-\bar{X})Y_i}{\Sigma(X_i-\bar{X})^2})\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X})Y_i)\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X})(\alpha + \beta X_i+\varepsilon_i))\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X}(\alpha + \beta X_i)+\Sigma(X_i-\bar{X})\varepsilon_i)\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}Var(\Sigma(X_i-\bar{X})\varepsilon_i)=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}\Sigma Var((X_i-\bar{X})\varepsilon_i)\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}(\Sigma (X_i-\bar{X})^2 Var(\varepsilon_i))\)

   \(=\frac{1}{(\Sigma(X_i-\bar{X})^2)^2}(\Sigma (X_i-\bar{X})^2 \sigma_\varepsilon^2)\)

   \(=\frac{\sigma_\varepsilon^2}{\Sigma(X_i-\bar{X})^2}\)

4. Derive the sampling variances of \(B\) in a simple regression.

   \(Var(A)=Var(\bar{Y}-B\bar{X})\)

   \(=Var(\bar{Y})+Var(B\bar{X})-2Cov(\bar{Y},B\bar{X})\)

   \(=Var(\bar{Y})+Var(B\bar{X})\)

   \(=Var(\frac{\Sigma Y_i}{n})+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}Var(\Sigma Y_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}\Sigma Var(Y_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}\Sigma Var(\alpha+\beta X_i+\varepsilon_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}\Sigma Var(\varepsilon_i)+\bar{X}^2Var(B)\)

   \(=\frac{1}{n^2}n\sigma_\varepsilon^2+\bar{X}(\frac{\sigma_\varepsilon^2}{\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{1}{n}+\frac{\bar{X}^2}{\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma(X_i-\bar{X})^2 +n\bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-2\bar{X}\Sigma X_i+\Sigma \bar{X}^2+n\bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-2\bar{X}\Sigma X_i+2n \bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-2\frac{\Sigma X_i}{n}\Sigma X_i+2n \bar{X}^2}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\sigma_\varepsilon^2(\frac{\Sigma X_i^2-\frac{2\Sigma X_i^2}{n}+2n\frac{\Sigma X_i^2}{n^2}}{n\Sigma(X_i-\bar{X})^2})\)

   \(=\frac{\sigma_\varepsilon^2 \Sigma X_i^2}{n\Sigma(X_i-\bar{X})^2}\)

c. Confidence Interval & Hypothesis Testing

As mentioned earlier, the fifth assumption, normality, allows us to apply statistical inference tests to assess the \(\beta\)s we estimate. Recall from the lecture that the residual standard error \(S_E\) (how closely the regression line we estimate fits the scatter of our data points) is:

\[\hat{\sigma } = SE(E_i) = \sqrt{\frac{\Sigma (E_i)^2}{n-2}} = \sqrt{\frac{\Sigma (Y_i-\hat{Y_i})^2}{n-2}}\]

And the standard error of the sample intercept (\(A\)) and slope (\(B\)) are:

\[SE(A) = \hat{\sigma }\sqrt{\frac{\Sigma X_i^2}{n\Sigma (X_i-\bar{X})^2}}\]

\[SE(B) = \frac{\hat{\sigma }}{\sqrt{\Sigma (X_i-\bar{X})^2}}\]

With standard errors, we can construct a \(100(1-\alpha)%\) confidence interval for our slope and perform hypothesis testing.

• Confidence Interval

  \(\beta=B+t_{\frac{\alpha}{2}}SE(B)\)

• Hypothesis Test

  Two-tailed test (by defualt, lm()function conducts two-tailedtest)

  \(H_0: \beta = 0\)

  \(H_1: \beta \neq 0\)

  \(t=\frac{b_1-\beta_1}{SE(b_1)}=\frac{b_1-0}{SE(b_1)}\)

  Compare this t-statistic with \(t_{\frac{\alpha}{2}, df=n-2}\) to see if it is larger or smaller than the critical value or not. If it is larger or smaller than the critical value, we can reject \(H_0\) in favor of \(H_1\). In addition to calculating t-statistics, we can also calculate p-value based the t-statistics and the critical values to perform hypothesis testing (p-value \(= 2*Pr(t \geq |t_c|)\) vs \(\alpha\)).

  One-tailed test (if your hypotheses are directional, you can also conduct one-tailed test)

  \(H_0: \beta = 0\)

  \(H_1: \beta < or > 0\)

  \(t=\frac{B-\beta_1}{SE(B)}=\frac{B-0}{SE(B)}\)

  When performing a one-tailed test, compare the t-statistic with \(t_{\alpha, df=n-2}\) without dividing \(\alpha\) by 2. P-value in one-tailed test is calculated as:

  p-value \(= Pr(t \geq or \leq t_c)\) vs \(\alpha\)

  Let’s run a simple example of it using the prestige data from Problem Set 1. The prestige data consists of 10 observations with 2 variables. The description of the variables are as follows:

  • prestige:The average prestige rating for the occupation.

  • education: The average number of years of education for occupational incumbents.

  # load data
  df <- read.csv("/Users/yu-shiuanhuang/Desktop/method-sequence/data/ps1_prestige.csv") 
  df

  ##    prestige education
  ## 1        82        86
  ## 2        83        76
  ## 3        90        92
  ## 4        76        90
  ## 5        90        86
  ## 6        87        84
  ## 7        93        93
  ## 8        90       100
  ## 9        52        87
  ## 10       88        86

  Let’s first calculate all the statistics we’ll need to perform hypothesis testing.

  1. Slope: \(B\)

     \(B=\frac{\Sigma(X_i-\bar{X})(Y_i-\bar{Y})}{\Sigma(X_i-\bar{X})^2}\)

  # calculate the slope b1
  B <- sum((df$education-mean(df$education))*(df$prestige-mean(df$prestige)))/
  sum((df$education-mean(df$education))^2)
  B

  ## [1] 0.3895028

  2. Intercept: \(A\)

     \(A=\bar{Y}-B\bar{X}\)

  # calculate the intercept b0
  A <- mean(df$prestige)-B*mean(df$education)
  A

  ## [1] 48.82376

  3. Residual Standard Error: \(\hat{\sigma} = SE(E_i)\)

     \(\hat{\sigma } = SE(E_i) = \sqrt{\frac{\Sigma (E_i)^2}{n-2}} = \sqrt{\frac{\Sigma (Y_i-\hat{Y_i})^2}{n-2}}\)

  # generate predicted Y
  df$pred <- A + B*df$education
  
  #calculate the residual standard error 
  se <- sqrt(sum((df$prestige-df$pred)^2)/(nrow(df)-2))
  se

  ## [1] 12.46986

  4. Standard Error of the Estimated Intercept and Slope: \(SE(A)\) & \(SE(B)\)

     \(SE(A) = \hat{\sigma }\sqrt{\frac{\Sigma X_i^2}{n\Sigma (X_i-\bar{X})^2}}\)

     \(SE(B) = \frac{\hat{\sigma }}{\sqrt{\Sigma (X_i-\bar{X})^2}}\)

  # calculate the standard error of the estimated intercept b0
  se_A <- se*sqrt(sum(df$education^2)/(nrow(df)*sum((df$education-mean(df$education))^2)))
  se_A

  ## [1] 57.80998

  # calculate the standard error of the estimated slope b1
  se_B <- se/sqrt(sum((df$education-mean(df$education))^2))
  se_B

  ## [1] 0.6554015

  Now, let’s perform two-tailed hypothesis testing of the slope.

  \(H_0: \beta = 0\)

  \(H_1: \beta \neq 0\)

  \(t=\frac{B-\beta_1}{SE(B)}=\frac{B-0}{SE(B)}=\frac{0.3895028-0}{0.6554015}=0.5942964\)

  # calculate the t-statistic
  t <- B/se_B
  t

  ## [1] 0.5942964

  What are the t critical values in both tails? Let’s say the level of significance, \(\alpha\), is 0.05 in this case.

  \(|t_{\frac{\alpha}{2}, df=n-2}|=|t_{\frac{0.05}{2}, df=10-2}|=2.306004\)

  # find t critical value
  qt(p = 0.025, df = 8, lower.tail = TRUE)

  ## [1] -2.306004

  qt(p = 0.025, df = 8, lower.tail = FALSE)

  ## [1] 2.306004

  Now, let’s compare our t-statistic of \(b_1\) with the t critical value.

  \(t=0.5942964 < |t_{\frac{0.05}{2}, df=8}|=2.306004\)

  We can also calculate the p-value.

  p-value \(= 2*Pr(t \geq 0.5942964)=0.5687352 > 0.05\)

  # find p-value
  2*pt(q = 0.5942964, df = 8, lower.tail = FALSE)

  ## [1] 0.5687352

  These all suggest that we cannot reject \(H_0\), which means that the positive effect of education on one’s prestige is not statistically significant.

  Let’s also try perform one-tailed hypothesis testing for this case.

  \(H_0: \beta = 0\)

  \(H_1: \beta > 0\)

  \(t=\frac{B-\beta_1}{SE(b_1)}=\frac{B-0}{SE(B)}=\frac{0.3895028-0}{0.6554015}=0.5942964\)

  What is the t critical values in the right tail?

  \(t_{\alpha, df=n-2}=t_{0.05, df=10-2}=1.859548\)

  # find t critical value
  qt(p = 0.05, df = 8, lower.tail = FALSE)

  ## [1] 1.859548

  Now again, let’s compare our t-statistic of \(B\) with the t critical value.

  \(t=0.5942964 < t_{0.05, df=8}=1.859548\)

  And the p-value in the one-tailed test is:

  p-value \(=Pr(t \geq 0.5942964)= 0.2843676 > 0.05\)

  # find p-value
  pt(q = 0.5942964, df = 8, lower.tail = FALSE)

  ## [1] 0.2843676

  In R, lm() does all the above calculation for us!

  # simple regression
  fit <- lm(prestige ~ education, data = df)
  summary(fit)

  ## 
  ## Call:
  ## lm(formula = prestige ~ education, data = df)
  ## 
  ## Residuals:
  ##      Min       1Q   Median       3Q      Max 
  ## -30.7105   0.3157   4.9580   5.6238   7.9525 
  ## 
  ## Coefficients:
  ##             Estimate Std. Error t value Pr(>|t|)
  ## (Intercept)  48.8238    57.8100   0.845    0.423
  ## education     0.3895     0.6554   0.594    0.569
  ## 
  ## Residual standard error: 12.47 on 8 degrees of freedom
  ## Multiple R-squared:  0.04228,    Adjusted R-squared:  -0.07743 
  ## F-statistic: 0.3532 on 1 and 8 DF,  p-value: 0.5687

d. Evaluating Least Squares Fit

We understand that there will be variability in the variable prestige across different occupations, with some having high prestige and others having low prestige. Some of this variability can be attributed to factors not included as variables in our regression model, such as economic conditions and other unmeasured variables. Additionally, differences in the variable education may explain some of the differences in prestige across occupations. For instance, occupations with a higher average number of years of education among incumbents tend to have higher prestige. However, it’s important to quantify how much of the variability in prestige can be explained by the variables included in our regression model, and how much remains unexplained by other factors that we have not accounted for.

1. Explained Variability: Regression Sum of Squares (RegSS)

   The Regression Sum of Squares (RegSS) (sometimes referred to as the Sum of Squares Explained, SSE) is a measure of the variability in the outcome variable that is explained by the explanatory variables, i.e. the x-variables in your regression. It is given by the following sum:

   \(RegSS = \sum_{i=1}^{n}(\hat{y_i}-\bar{y})^2\)

2. Residual or Unexplained Variability: Residual Sum of Squares (RSS)

   The Residual Sum of Squares (RSS) is a measure of the variability in the outcome variable that is not explained by your regression, and therefore is due to all the other factors that affect prestige besides education. It is given by the following sum:

   \(RSS = \sum_{i=1}^{n}(y_i-\hat{y_i})^2 = \sum_{i=1}^{n}e_i^2\)

3. Total Variability: Total Sum of Squares (TSS)

   You can show mathematically that \(RSS + RegSS\) is equal to the following expression, which is referred to as the Total Sum of Squares (TSS):

   \(TSS = \sum_{i=1}^{n}(y_i-\bar{y_i})^2 = \sum_{i=1}^{n}(y_i-\hat{y_i})^2 + \sum_{i=1}^{n}(\hat{y_i}-\bar{y})^2 = RSS + RegSS\)

4. Coefficient of Determination: R-Squared value

   The coefficient of determination, sometimes referred to as the R-Squared value, is a measure of what percentage of the variability in your outcome variable is explained by your explanatory variables. It is given by the expression,

   \(R^2 = \frac{RegSS}{TSS}\)

Discussion II

For today’s discussion, we will review concepts regarding how to find a set of \(\hat{\beta}\)s that can minimize the squared residual errors of our multivariate linear model using matrix algebra.

1. Multivariate Linear Regression
2. Find Coefficients using Matrix in R

1. Multivariate Linear Regression

The central difference between simple and multiple linear regression is that the slope coefficients for the explanatory variables in multiple regression are partial coefficients. That is, it represents the “effect” on the response variable of a one-unit increment in the corresponding explanatory variable, holding constant the value of the other explanatory variable.

\[Y_i=\hat{\beta_0}+\hat{\beta_1}X_{i1}+\hat{\beta_2}X_{i2}+...+\hat{\beta_k}X_{ik}+\varepsilon_i\]

Recall that the goal of the OLS approach is to minimize the discrepancy between our observed \(Y_i\) and estimated \(\hat{Y_i}\). To find a set of estimated \(\beta\)s that can minimize this discrepancy is an optimization problem. All we have to do is to take the partial derivatives of the above equation with respect to each coefficient, set them equal to zero, and then solve it. However, as shown in the lecture, adding more and more \(X\) variables to our linear regression model makes the calculation more and more tedious. In order to find the estimated \(\beta\)s more efficiently, we use matrix algebra to help us do the calculation.

To do this, let’s first think of our data set as a matrix. The below equations are how each observation in our data set generated.

\[Y_1=\hat{\beta_0}+\hat{\beta_1}X_{11}+\hat{\beta_2}X_{12}+...+\hat{\beta_k}X_{1k}+\varepsilon_1\]

\[Y_2=\hat{\beta_0}+\hat{\beta_1}X_{21}+\hat{\beta_2}X_{22}+...+\hat{\beta_k}X_{2k}+\varepsilon_2\]

.

.

.

\[Y_n=\hat{\beta_0}+\hat{\beta_1}X_{n1}+\hat{\beta_2}X_{n2}+...+\hat{\beta_k}X_{nk}+\varepsilon_n\]

Rewrite them into a matrix form.

The system of equations is summarized as the y vector equal to the product of the X matrix of data times the \(\hat{\beta}\) vector of parameters plus the vector of disturbances, \(\varepsilon\). That said, the fitted linear model we are trying to find is: \(y = X\hat{\beta}+\varepsilon\). In the lecture, Lauren showed how to apply what we did when finding a pair of intercept and slope in simple linear regression to solve the optimization problem(i.e., minimize RSS) to find the coefficient vector \(\hat{\beta}\) using matrix algebra.

\[ \begin{aligned} S(\hat{\beta}) &= \sum E_i^2 = \varepsilon\varepsilon' = (y- X\hat{\beta})'(y- X\hat{\beta})\\ &= y'y-y'X\hat{\beta}-\hat{\beta}'X'y+\hat{\beta}'X'X\hat{\beta} \\ &= y'y-(2y'X)\hat{\beta}+\hat{\beta}'(X'X)\hat{\beta} \end{aligned} \]

Instead of taking derivative with respect to each coefficient one at a time, we can take derivative directly with respect to the coefficient vector \(\hat{\beta}\).

\[ \begin{aligned} S(\hat{\beta}) &= \sum E_i^2 = \varepsilon\varepsilon' = (y- X\hat{\beta})(y- X\hat{\beta})'\\ &= y'y-y'X\hat{\beta}-\hat{\beta}'X'y+\hat{\beta}'X'X\hat{\beta} \\ &= y'y-(2y'X)\hat{\beta}+\hat{\beta}'(X'X)\hat{\beta} \end{aligned} \] \[ \begin{aligned} \partial S(\hat{\beta}) &= \frac{\partial}{\partial \hat{\beta}}(y'y-(2y'X)\hat{\beta}+\hat{\beta}'(X'X)\hat{\beta}) \\ 0 &=0-2Xy'+ 2\hat{\beta}X'X \\ 2\hat{\beta}X'X &=2Xy'\\ \hat{\beta} &= \frac{Xy'}{X'X} \\ \hat{\beta} &= (X'X)^{-1}Xy' \end{aligned} \] Therefore, the least squares estimator is:

\[\hat{\beta}=(X'X)^{-1}X'y\]

2. Find Coefficients using Matrix in R

Let’s use the Anscombe.txt data set as an example to find the estimated coefficients by using the above matrix equation.

• education = per capita education expenditures, dollars
• income = per capita income, dollars
• under18 = proportion under 18 years old, per 1000
• urban = proportion urban, per 1000

# load data
df2 <- read.delim("/Users/yu-shiuanhuang/Desktop/method-sequence/data/Anscombe.txt", sep = "")
head(df2)

##    education income under18 urban
## ME       189   2824   350.7   508
## NH       169   3259   345.9   564
## VT       230   3072   348.5   322
## MA       168   3835   335.3   846
## RI       180   3549   327.1   871
## CT       193   4256   341.0   774

# convert each variable to a vector
edu <- df2$education
intercept <- rep(1, nrow(df2)) # remember to crease this vector with 1s!
inc <- df2$income
un18 <- df2$under18
urb <- df2$urban
X <- cbind(intercept, inc, un18, urb)
class(X) # this is a matrix

## [1] "matrix" "array"

We now have a vector edu and a matrix \(X\) containing our intercepts and our three explanatory variables (income, under18, urban). Let’s estimate our best fit parameters \(\beta_0, \beta_1,\), \(\beta_2\), and \(\beta_3\) by applying the estimator we derived in the lecture: \(\hat{\beta} = (X'X)^{-1}X'y\).

beta.hat <- solve(t(X) %*% X) %*% t(X) %*% edu

t(beta.hat)

##      intercept        inc      un18        urb
## [1,] -286.8388 0.08065325 0.8173377 -0.1058062

Compare with the canned routine in R to confirm the matrix calculation was correct.

mod <- lm(edu ~ X-1) # remember to remove the first column in the X matrix bc that's the column with all the 1s we created just for the algebra calculation

beta.canned <- mod$coefficients
beta.canned

##    Xintercept          Xinc         Xun18          Xurb 
## -286.83876273    0.08065325    0.81733774   -0.10580623

How to come up with the predicted \(\hat{y}\)? \[\hat{y}=Xb\]

eduhat <- beta.hat[1] + df2$income*beta.hat[2] + 
    df2$under18*beta.hat[3] + df2$urban*beta.hat[4]
eduhat

##  [1] 173.8168 199.0526 211.7006 207.0077 174.5936 253.2396 223.9983 210.5846
##  [9] 179.8788 206.2476 213.3513 231.5918 233.2380 209.4856 211.0236 196.9737
## [17] 176.3859 188.1518 199.2828 195.3128 187.4341 250.0855 231.5108 252.0303
## [25] 182.3619 139.8510 169.8280 162.6433 176.5621 159.9772 156.6485 139.1353
## [33] 135.8968 148.7554 135.1561 174.0988 143.0524 175.0569 195.4585 171.6380
## [41] 205.2510 192.1739 197.6281 189.7953 190.1857 261.4657 212.7525 181.6204
## [49] 221.7760 355.7228 221.5298

How to come up with the residuals? \[e=y-\hat{y}=y-Xb\]

e <- edu - eduhat
e

##  [1]  15.1831901 -30.0526055  18.2993676 -39.0077435   5.4064124 -60.2396370
##  [7]  37.0016640   3.4153951  21.1211678 -34.2476465 -19.3513498 -42.5918115
## [13]  -0.2380450  -0.4855509  50.9763644  37.0263340   0.6140641 -11.1518234
## [19] -12.2828447 -47.3127738   8.5658819  -2.0854980  15.4891839  -6.0302784
## [25]  -2.3619147   9.1490037 -14.8279988 -13.6433457 -20.5621494  31.0227604
## [31] -16.6485112  -2.1352962 -23.8967859 -18.7553867  -1.1561315 -12.0987770
## [37]  -8.0523620 -20.0569429  42.5415211  -1.6380384  32.7489708  -0.1738631
## [43]  29.3718741  17.2047442  10.8143055 -36.4656906   2.2475110  51.3796304
## [49]  51.2240417  16.2771794  -9.5297657

Now calculate the standard errors of the vector \(\hat{\beta}\). Recall the error variance is unknown but we have available the unbiased estimator \(S^2_E = (e'e)/(n-k-1)\), based on the residuals of our model fit. We can extract these residuals, calculate the estimator, and then use it to calculate the variance of the least squares coefficients \(\hat{V(b)} = S^2_E (X'X)^{-1}\).

• As a quick refresher of concepts: the variance is a measure of a random variable’s “spread” or variation around its mean (a.k.a. its expected value), while the co-variance measures how correlated are the variations of two random variables with each other.

• The variance-covariance matrix is a square matrix (i.e. it has the same number of rows and columns). The elements of the matrix that lie along its main diagonal (i.e. the one that goes from top-left to bottom-right contain the variances while all other elements contain the co-variances). Thus, the variance-covariance matrix of the fitted coefficients of a regression model contains the variances of the fitted model’s coefficient estimates along its main diagonal, and it also contains the pair-wise co-variances between coefficient estimates in the non-diagonal elements.

• Why do we not only calculate the variance of regression coefficients but also the convariance of them?

  The covariance matrix can be very useful for certain model diagnostics. If two variables are highly correlated, one way to think about it is that the model is having trouble figuring out which variable is responsible for an effect (because they are so closely related). This can be helpful for a whole variety of cases, such as choosing subsets of covariates to use in a predictive model; if two variables are highly correlated, you may only want to use one of the two in your predictive model.

S2E <- as.numeric((t(e)%*%e)/(nrow(df2)-3-1))
S2E

## [1] 712.5394

v.beta.hat <- S2E*solve(t(X)%*%X)
v.beta.hat # the diagonal numbers are the variance of each coefficient

##              intercept           inc          un18           urb
## intercept 4214.5975090 -1.849526e-01 -9.7493677311 -0.1582902773
## inc         -0.1849526  8.646281e-05  0.0001387409 -0.0002162611
## un18        -9.7493677  1.387409e-04  0.0255327479  0.0002084764
## urb         -0.1582903 -2.162611e-04  0.0002084764  0.0011752674

sqrt(diag(v.beta.hat)) 

##    intercept          inc         un18          urb 
## 64.919931523  0.009298538  0.159789699  0.034282173

Compare with the canned routine in R to see if we got the same standard errors of coefficients.

	OLS
Variables	Estimates	S.E.	C.I. (95%)
Xintercept	-286.84 ***	64.92	-417.44 – -156.24
Xinc	0.08 ***	0.01	0.06 – 0.10
Xun18	0.82 ***	0.16	0.50 – 1.14
Xurb	-0.11 **	0.03	-0.17 – -0.04
Observations	51
R2 / R2 adjusted	0.984 / 0.982
p<0.05   ** p<0.01   *** p<0.001

With the standard errors of each coefficient, we can perform hypothesis testing and construct confidence intervals to determine the significance and reliability of the estimated effects (similar to what we did in the case of simple linear regression).

• What would happen when there is perfect multicollinearity in our data (i.e. when two or more independent variables are perfectly correlated)?

  # convert each variable to a vector
  edu <- df2$education
  intercept <- rep(1, nrow(df2)) # remember to crease this vector with 1s!
  inc <- df2$income
  un18 <- df2$under18
  urb <- df2$urban
  pmc <- 2*inc
  cor(pmc, inc) # pmc is perfectly correlated with inc

  ## [1] 1

  X <- cbind(intercept, inc, un18, urb, pmc)
  class(X) # this is a matrix

  ## [1] "matrix" "array"

  beta.hat <- solve(t(X) %*% X) %*% t(X) %*% edu

  ## Error in solve.default(t(X) %*% X): Lapack routine dgesv: system is exactly singular: U[5,5] = 0

  As you can see in the result of the above chunk, you cannot calculate the matrix algebra when there is perfect multicollinearity. However, when you run lm() in R, R will automatically delete the variable that is perfectly correlated with the other variable(s) in the dataset.

  mod2 <- lm(edu ~ X-1) # remember to remove the first column in the X matrix bc that's the column with all the 1s we created just for the algebra calculation

  	OLS
  Variables	Estimates	S.E.	C.I. (95%)
  Xintercept	-286.84 ***	64.92	-417.44 – -156.24
  Xinc	0.08 ***	0.01	0.06 – 0.10
  Xun18	0.82 ***	0.16	0.50 – 1.14
  Xurb	-0.11 **	0.03	-0.17 – -0.04
  Observations	51
  R2 / R2 adjusted	0.984 / 0.982
  p<0.05   ** p<0.01   *** p<0.001

  If there is high multicollinearity between variables but not perfectly correlated, the matrix algebra still works, and lm() will still spit out the estimates of the highly correlated variable. However, the standard errors of coefficients will be very high.

  # convert each variable to a vector
  edu <- df2$education
  intercept <- rep(1, nrow(df2)) # remember to crease this vector with 1s!
  inc <- df2$income
  un18 <- df2$under18
  urb <- df2$urban
  set.seed(1234)
  pmc <- 2*inc + rnorm(nrow(df2), 5, 1)
  cor(pmc, inc) # pmc is highly but not perfectly correlated with inc

  ## [1] 0.9999997

  X <- cbind(intercept, inc, un18, urb, pmc)
  class(X) # this is a matrix

  ## [1] "matrix" "array"

  beta.hat <- solve(t(X) %*% X) %*% t(X) %*% edu
  
  t(beta.hat)

  ##      intercept     inc      un18        urb       pmc
  ## [1,] -290.6956 -1.4651 0.8188972 -0.1038287 0.7726427

  eduhat <- beta.hat[1] + df2$income*beta.hat[2] + 
  df2$under18*beta.hat[3] + df2$urban*beta.hat[4] + pmc*beta.hat[5]
  e <- edu - eduhat
  S2E <- as.numeric((t(e)%*%e)/(nrow(df2)-4-1))
  v.beta.hat <- S2E*solve(t(X)%*%X)
  sqrt(diag(v.beta.hat))

  ##   intercept         inc        un18         urb         pmc 
  ## 69.35278085  9.01847688  0.16172179  0.03651251  4.50787061

  	OLS
  Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
  Xintercept	-286.84 ***	64.92	-417.44 – -156.24	-290.70 ***	69.35	-430.30 – -151.10
  Xinc	0.08 ***	0.01	0.06 – 0.10	-1.47	9.02	-19.62 – 16.69
  Xun18	0.82 ***	0.16	0.50 – 1.14	0.82 ***	0.16	0.49 – 1.14
  Xurb	-0.11 **	0.03	-0.17 – -0.04	-0.10 **	0.04	-0.18 – -0.03
  Xpmc				0.77	4.51	-8.30 – 9.85
  Observations	51			51
  R2 / R2 adjusted	0.984 / 0.982			0.984 / 0.982
  p<0.05   ** p<0.01   *** p<0.001

Discussion III

For today’s discussion, we will review some diagnostic techniques to detect influential data points that could affect our model and ensure our model adheres to Gauss-Markov assumptions, addressing any violations as needed.

1. Unusual and Influential Data
2. Diagnostics for Errors and Linearity

1. Unusual and Influential Data

When dealing with observational or experimental data, it’s common to encounter outliers, which are data points that behave differently than the rest of the data for some reason. Outliers can greatly influence the analysis and interpretation of the data. In today’s discussion, we’ll examine what outliers are, how they can affect data analysis and interpretation, and use simple linear regression as an example to explore situations where individual data points significantly impact the model. Lastly, we will also

1.1 Leverage, discrepancy and influence

Although many researchers discuss outliers in a broad sense, it’s important to recognize that there are different types of unusual data points. A data point may be unusual in terms of its predictor behavior (the \(x\)-value in simple regression), its outcome (the \(y\)-value in simple regression), or both.

• leverage: A data point that differs significantly from the rest of the dataset in terms of its predictor behavior is known as a leverage point. In the case of simple linear regression, this means that the x-value of the data point is much higher or lower than the mean of the predictor variable.

• discrepancy: When a data point has an uncommon y-value for its corresponding x-value, it is said to have high discrepancy. Within the context of regression analysis, such an unusual point is called an outlier.

• influence: Having either a high leverage or a high discrepancy does not necessarily make a data point influential in a linear model. In fact, a single data point’s influence is determined by leverage × discrepancy. Therefore, having only high leverage or high discrepancy may not be enough to significantly alter the model parameters.

The figure below illustrates the three characteristics discussed above. In each of the three panels, the red line represents the line of best fit without the data point in question (marked by the triangle), while the blue line shows the line of best fit with it.

• In panel A, the data point with the triangle has a high leverage - its \(x\)-value is much higher than the rest.

• In panel B, it has a high discrepancy - it lies far away from the line of best fit.

  However, neither of the above characteristics alone exerts a lot of influence on our model parameters, as the red line does not diverge significantly from the blue line.

• In panel C, we see a data point that has both high leverage and high discrepancy, and as a result, it exerts high influence: the blue line is noticeably different from the red line.

Example Data

For today’s discussion, we will use a data set based on an example in Field, Miles, and Field (2012). The dataset includes eight samples, where the \(x\)-value represents the number of pubs within a borough (district) of London, and the \(y\)-value represents the number of deaths in that borough over a specific period of time. We are interested in investigating the relationship between the number of deaths and the number of pubs in each borough.

pubs <- data.frame(name = c(1, 2, 3, 4, 5, 6, 7, 8), 
                   pubs = c(10, 25, 40, 55, 70, 85, 100, 500), 
                   deaths = c(1043.822, 2086.934, 2951.086, 3992.459, 5088.003, 
                              6095.645, 6923.497, 10000.000))
pubs

##   name pubs    deaths
## 1    1   10  1043.822
## 2    2   25  2086.934
## 3    3   40  2951.086
## 4    4   55  3992.459
## 5    5   70  5088.003
## 6    6   85  6095.645
## 7    7  100  6923.497
## 8    8  500 10000.000

Let’s plot the line of best fit of this data:

library(tidyverse)
ggplot(data = pubs, aes(x = pubs, y = deaths, label = name)) +
    geom_point(alpha = 0.8) +
    geom_text(hjust = 2, vjust = 0) +
    geom_smooth(method = lm, se = FALSE, colour = "darkorange") +
    xlab("Number of Pubs") + 
    ylab("Deaths") +
    theme_bw()

Both the raw data and the linear regression plot indicate that data point 8 is significantly different from the other samples. This particular data point not only has a much higher number of pubs compared to the others, but it also appears to have a distinct relationship between the number of pubs and the number of deaths, unlike the other boroughs.

Since we have some concerns about this specific data point, let’s determine whether it exhibits a) high leverage, b) high discrepancy, and c) high influence.

1.2 Assessing Leverage

Recall that leverage quantifies the degree to which a predictor value differs from the other predictor values. In the case of simple linear regression, we can compute the distance between each predictor point (\(X_i\)) and the mean predictor value (\(\bar{X}\)) to measure leverage.

A standardized version of this distance is called hat-value and denoted by \(h_i\) (for \(i = {1,..., n}\)):

\[h_i = \frac{1}{n}+\frac{(X_i-\bar{X})^2}{\sum_{j=1}^{n}(X_j-\bar{X})^2}\]

The average hat value (\(\bar{h}\)) is defined as \(\frac{k+1}{n}\), in which \(k\) is the number of predictors and \(n\) the number of observations. Values of \(h\) are bound between \(\frac{1}{n}\) and 1, with 1 denoting highest leverage (highest distance from mean).

Based on our sample data, it’s obvious that data point 8 has the highest leverage among all the points. Let’s apply the formula above to compute the leverage value (\(h_8\)) for this data point.

# number of cases
n = nrow(pubs)

# distance to mean of point 8
numerator = (pubs$pubs[8] - mean(pubs$pubs))^2

# distance to mean of all the other points
denominator = sum((pubs$pubs - mean(pubs$pubs))^2)

# putting it together
h_8 = 1/n + numerator/denominator

h_8

## [1] 0.969302

The resulting hat-value is 0.969302. This is quite high - in fact, it’s very close to 1, the highest possible value.

You don’t need to manually calculate all of the hat values because R offers a convenient hatvalues function that can be applied to any linear model. To take advantage of this feature, we will begin by fitting a simple linear model using lm and then extracting the hat value of each data point.

# fitting linear model
mod <- lm(deaths ~ pubs, data = pubs)

# getting hatvalues and plotting them
plot(hatvalues(mod))

Remember that leverage alone does not mean a point exerts high influence, but it certainly means it’s worth investigating. Hat values are open to interpretation, but a cut-off value that is common is twice or three times the average \(\bar{h}\), meaning anything above that value should be looked at closer.

Let’s add a cut-off line into the above hatvalue figure.

# the average of hat-values
mean_h <- (1+1)/8

# adding a cut-off line
plot(hatvalues(mod))
abline(h = 2*mean_h, lty = 2)
abline(h = 3*mean_h, lty = 2)

As illustrated in the figure above, we can see that the value of \(h_8\) is greater than twice and three times the mean hat values (\(\bar{h}\)). Based on this observation, we can safely conclude that the value of \(h_8\) is highly unusual. In this case, \(h_8\) is definitely unusual!

After having assessed leverage, let’s look at discrepancy. How unusual is the \(y\)-value given its \(x\)-value?

1.3 Assessing Discrepancy

As we learned earlier, data points that do not fit well to the linear regression line are referred to as outliers, or high-discrepancy points. Typically, we assess the fit of the regression line using residuals, which measure the distance between a predicted value and the actual value. While we already know that point 8 has high leverage, it is worth noting that the line of best fit is quite close to its predicted value (see the left panel). Additionally, when we examine the residuals, point 8 does not appear to be further from the regression line than any other points (see the right panel).

Instead, we can look at standardized or studentized residuals.

• Standardized residuals are formed by calculating:

  \[E_i'=\frac{E_i}{S_E\sqrt{1-h_i}}\]

  where \(S_E=\sqrt{\frac{\sum E_i^2}{n-k-1}}\)

  Let’s calculate the standardized residual for our example data by hand (the 8th data point is the one we’re interested in):

  # residual for data point in original model
  Ei <- as.numeric(residuals(mod)[8])
  
  # estimate of sigma (standard deviation) for residuals
  S_E <- summary(mod)$sigma
  
  # hatvalue for point 8
  hi <- as.numeric(hatvalues(mod)[8])
  
  # putting it together
  Eprime <- Ei/(S_E*sqrt(1-hi))
  Eprime

  ## [1] -2.447433

  You can also just let R do the maths by calling rstandard on the original model:

  rstandard(mod)

  ##           1           2           3           4           5           6 
  ## -1.44154418 -0.89769067 -0.48020238  0.04464274  0.59853447  1.09293034 
  ##           7           8 
  ##  1.47205817 -2.44743280

  Note that the last value is the same as the one we calculated by hand! Let’s compare the residuals and the studentized residuals in plots side by side:

  

  The standardized residuals provide a better measure of discrepancy and can reveal outliers that were not detected by the normal residual plot. In our case, by looking at the standardized residual plot below, we can confirm that point 8 has a high discrepancy compared to the other points.

• Studentized residuals are calculated by fitting a model without the case for which the residual is calculated, and then scaling the resulting residual (\(E_i\)) by an estimate of the standard deviation of the residuals (\(S_{E(-i)}\)) and the point’s hat value (\(h_i\)):

  \[E_i^*=\frac{E_i}{S_{E(-i)}\sqrt{1-h_i}}\]

  Let’s calculate the studentized residual for our example data by hand (the 8th data point is the one we’re interested in):

  # model excluding point 8
  mod.red <- lm(deaths ~ pubs, data = pubs[-8,]) 
  
  # residual for data point in original model
  Ei <- as.numeric(residuals(mod)[8])
  
  # estimate of sigma (standard deviation) for residuals
  S_E <- summary(mod.red)$sigma
  
  # hatvalue for point 8
  hi <- as.numeric(hatvalues(mod)[8])
  
  # putting it together
  Estar <- Ei/(S_E*sqrt(1-hi))
  Estar

  ## [1] -54.5284

  You can also just let R do the maths by calling rstudent on the original model:

  rstudent(mod)

  ##            1            2            3            4            5            6 
  ##  -1.62765333  -0.88075356  -0.44703731   0.04075983   0.56346508   1.11482846 
  ##            7            8 
  ##   1.68127219 -54.52839824

  Note that the last value is the same as the one we calculated by hand! Let’s compare the residuals and the studentized residuals in plots side by side:

  

  The studentized residuals reveal clearly that point 8 has a high discrepancy, while this was not possible to see from the normal residual plots. You can also perform Bonferroni adjustment to test the statistical significance of the \(E_8^*\) to check if we can confidently say that point 8 is an outlier. Observations are considered outliers if their Bonferroni p is less than .05.

  library(car)
  outlierTest(mod)

  ##   rstudent unadjusted p-value Bonferroni p
  ## 8 -54.5284         3.9231e-08   3.1385e-07

  According to the above test, we can say that point 8 can be considered an outlier (high discrepancy).

1.4 Assessing Influence

As mentioned earlier, having high leverage or discrepancy alone does not necessarily mean a data point will have a high influence on the linear regression model. However, since our point of interest (point 8) has both high leverage and discrepancy, it is likely to have high influence on the model.

We’ll use something called Cook’s D, which relies on the standardized residuals to assess influence. A similar measure is DFFITS, which instead is based on studentized residuals. Another measure we won’t go into are DFBETAS, which measures the influence on each individual parameter, instead of the overall model.

Observations that combine high leverage with large studentized residual exert substantial influence on regression coeffcients. Cook’s D statistic provides a summary index of influence.

\(D_i=\frac{E_i'^2}{k+1}\times \frac{h_i}{1-h_i}\)

where the first term, the standardized residual, is a measure of discrepancy and the second is a measure of leverage.The higher the leverage and residuals, the higher the Cook’s distance.

Several interpretations for Cook’s distance exist. There isn’t a universally accepted rule for cut off points.

• One interpretation is to investigate any point over \(\frac{4}{n}\), where \(n\) is the number of observations.
• Other scholars suggest that any “large” \(D_i\) should be investigated. How large is “too large”? The consensus seems to be that a \(D_i\) value of more than 1 indicates an influential value, but you may want to look at values above 0.5.

Let’s use influencePlot() in the car package to check Cook’s D statistics!

##      StudRes       Hat      CookD
## 1  -1.627653 0.1813863  0.2302244
## 7   1.681272 0.1256287  0.1556728
## 8 -54.528398 0.9693020 94.5671657

1.5 Added-variable Plot (AV Plot)

All the above techniques can also be applied to multivariate linear regression models. Additionally, we introduce another very useful influence graph for visually detecting outliers in multivariate linear models, known as an added-variable plot or partial-regression plot.

A limitation of using the typical residual plots in a multivariate linear regression model, such as plotting residuals against the values of a predictor variable, is that they may not fully illustrate the “additional contribution” of a predictor variable when considering other variables already in the model. Added-variable plots offer a more refined visualization, providing graphic insights into the marginal importance of a predictor variable, taking into account the presence of other variables in the model. The following materials are sourced from Fox’s slides and lecture.

• Let \(y_i^{(1)}\) represent the residuals from the least-squares regression of \(y\) on all of the \(x\)s with the exception \(x_1\), that is, the residuals from the fitted model:

  \[y_i = b_0^{(1)} + b_2^{(1)}x_{i2} + ...+ b_k^{(1)}x_{ik} + y_i^{(1)}\]

• Likewise, \(x_i^{(1)}\) are residuals from the least-squares regression of \(x_1\) on the other \(x\)s:

  \[x_{1i} = c_0^{(1)} + c_2^{(1)}x_{i2} + ...+ c_k^{(1)}x_{ik} + x_i^{(1)}\]

• The notation emphasizes the interpretation of the residuals \(y_i^{(1)}\) and \(x_i^{(1)}\) as the parts of \(y\) and \(x_1\) that remain when the linear dependence of these variables on \(x_2, ..., x_k\) is removed.

The AV plot for \(x_1\) is then the scatterplot of \(y_i^{(1)}\) versus \(x_i^{(1)}\), and we repeat the procedure for each \(x_j\), \(j =0, 1, ..., k\) (where \(x_0 = 1\)). In effect, the \((k+1)\)-dimensional scatterplot for \(y\) and \(x_1, ..., x_k\) is reduced to a sequence of \(k+1\) 2D AV plots. By doing so, the AV plots visualize leverage and influence on each of the regression coefficients.

The added-variable plot for \(x_1\) has the following very interesting properties:

• The slope of the least-squares simple-regression line of of \(y_i^{(1)}\) on \(x_i^{(1)}\) is the same as the least-squares slope \(b_1\) for \(x_1\) in the full multiple regression.

• The residuals from this simple regression are the same as the residuals \(e_i\) from the full regression.

• Consequently, the standard deviation of the residuals in the added-variable plot is the same as the multiple regression.

• The standard error of \(b_1\) from this simple regression is also the same as the standard error from the full regression

• Because the \(x_i^{(1)}\) are residuals, they are less variable than \(x_1\) if \(x_1\) is correlated with the other \(x\)s. The added-variable plot therefore show how collinearity can degrade the precision of the estimation by decreasing the conditional variation of an \(x\).

NOTE: You may find these concepts very familiar as in POL 212, we have actually covered the Frisch-Waugh-Lovell Theorem (you can review the materials of Discussion IV in POL 212 Discussion), in which we showed how to decompose a regression of \(Y\) on a set of variables \(X\) into two pieces. The coefficient of one of the variables in multiple linear regression model can be obtained by netting off effect of other variables in the regression model from both dependent and independent variables. The FWL Theorem shows that it is important to partial one set of covariates out of the other in order to obtain the ceteris paribus effect. For example, in order to get the effect of \(X1\) on \(Y\), \(X2\) needs to be held constant, which involves removing the portion of \(X1\) that is correlated with \(X2\). Likewise, when identifying outliers among explanatory variables in a multivariate linear regression model, it’s crucial to first hold all other variables constant.

Example

For this example, we will use the Duncan data frame, which has 45 rows and 4 columns. This is a dataset on the prestige and other characteristics of 45 U.S. occupations in 1950.

library(car)
library(carData)

knitr::kable(head(Duncan)) # from carData package

	type	income	education	prestige
accountant	prof	62	86	82
pilot	prof	72	76	83
architect	prof	75	92	90
author	prof	55	90	76
chemist	prof	64	86	90
minister	prof	21	84	87

These are the description of each variable:

• type: Type of occupation. A factor with the following levels: prof, professional and managerial; wc, white-collar; bc, blue-collar.

• income: Percentage of occupational incumbents in the 1950 US Census who earned $3,500 or more per year (about $36,000 in 2017 US dollars).

• education: Percentage of occupational incumbents in 1950 who were high school graduates (which, were we cynical, we would say is roughly equivalent to a PhD in 2017).

• prestige: Percentage of respondents in a social survey who rated the occupation as “good” or better in prestige.

Let’s begin by regressing prestige on all explanatory variables income and education, then plot the added-variable plots to visually inspect for any potential outliers.

m.duncan <- lm(prestige ~ income + education, data = Duncan)

	OLS
Variables	Estimates	S.E.	C.I. (95%)
(Intercept)	-6.06	4.27	-14.69 – 2.56
income	0.60 ***	0.12	0.36 – 0.84
education	0.55 ***	0.10	0.35 – 0.74
Observations	45
R2 / R2 adjusted	0.828 / 0.820
p<0.05   ** p<0.01   *** p<0.001

avPlots(m.duncan)

Note that the slope of the blue regression line in both plots corresponds to the estimated coefficient of income and education in our multivariate regression model.

In the AV plots, we can observe the impact (or lack thereof) of these unusual cases on the regression coefficients. Let’s first examine the left-hand side of the figure. We notice that the case of the minister pulls up on the left side, while the case of the conductor pulls down on the right side. This influential pair is widely separated, affecting the coefficient of income by making it smaller than it would otherwise be. The case of railroad engineers appears as a very high leverage point due to their unusual income given their education level, but it aligns more or less with the trend observed throughout the data. Similarly, the case of the reporter also follows the overall trend.

Similarly, in the education AV plots, we observe a similar trend except for the effect of the minister’s case, which pulls up on the left side, leading to a larger coefficient of education than it would otherwise be. Again, the cases of railroad engineers and reporters seem to align with the general trend of the data, requiring less concern.

Overall, based on the AV plots, we may speculate that the case of the minister could be an influential data point affecting our model’s estimation. We can proceed by updating our original regression model and then comparing the coefficients between the models to validate this hypothesis.

m.duncan2 <- lm(prestige ~ income + education, 
                data = subset(Duncan, !(rownames(Duncan) %in% c("minister"))))

	Model 1			Model 2
Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
(Intercept)	-6.06	4.27	-14.69 – 2.56	-6.63	3.89	-14.48 – 1.22
income	0.60 ***	0.12	0.36 – 0.84	0.73 ***	0.12	0.50 – 0.97
education	0.55 ***	0.10	0.35 – 0.74	0.43 ***	0.10	0.24 – 0.63
Observations	45			44
R2 / R2 adjusted	0.828 / 0.820			0.856 / 0.849
p<0.05   ** p<0.01   *** p<0.001

Indeed, after excluding the case of the minister, the coefficient of income increases from 0.6 to 0.73, while the coefficient of education decreases from 0.55 to 0.43.

We can also run the influencePlot() to check if the Cook’s distance of the minister’s case is indeed very problematic or not.

##                StudRes        Hat      CookD
## minister     3.1345186 0.17305816 0.56637974
## reporter    -2.3970224 0.05439356 0.09898456
## conductor   -1.7040324 0.19454165 0.22364122
## RR.engineer  0.8089221 0.26908963 0.08096807

2. Diagnostics for Linearity and Errors

2.1 Violation of Linearity Assumption

Recall the linearity assumption is that the expected (mean) value of the disturbance term is 0 (i.e., \(E(u_i)=0\)). The average regression error is 0 everywhere implies that the regression surface captures the dependency of the conditional mean of \(y\) on \(x\)s. Violating the linearity assumption implies “that the model fails to capture the systematic pattern of relationship between the response and explanatory variables” (Fox, 2016: 307). For example, a partial relationship specified to be linear may be nonlinear, or two explanatory variables specified to have additive partial effects may interact in determining \(y\). While a fitted model may sometimes provide a useful approximation to the true regression surface \(E(y)\), there are instances where the model can be highly misleading. You can think of nonlinearity (fitting the wrong equation to the data) as potentially the most serious problem with a regression model.

When performing a linear regression with a single explanatory variable, a scatterplot of the response variable against the explanatory variable provides a good indication of the nature of the relationship. If there is more than one explanatory variable, things become more complicated. Although it can still be useful to generate scatterplots of the response variable against each of the explanatory variables, this does not take into account the effect of the other explanatory variables in the model. After all, our interest centers on the partial relationship between \(y\) and each \(x\), controlling for the other \(x\)s, not on the marginal relationship between y and a single \(x\).

Furthermore, while plotting residuals against each \(x\) is helpful for detecting departures from linearity, residual plots cannot distinguish between monotone and nonmonotone nonlinearity.

In the above figure from the Fox textbook, case (a) might be modeled by \(y = \beta_0 + \beta_1\sqrt{x} + \varepsilon\) (a transformation of \(x\)), but case (b) cannot be linearied by a transformation of \(x\) and might instead be dealt with a quadratic regression, \(y = \beta_0 + \beta_1x + \beta_2x^2 + \varepsilon\).

Component + residual (C+R) plots, an extension of partial residual plots, are a good way to see if each predictor you include in the multivariate linear model have a linear relationship to the dependent variable. A partial residual plot essentially attempts to model the residuals of one predictor against the dependent variable. The partial residual for the \(j\)th explanatory variables is defined as:

\[E_{i}^{(j)}= E_i+B_jX_{ij}\]

In words, add back the linear component of the partial relationship between \(Y\) and \(X_j\) to the least-squares residuals, which may include an unmodeled nonlinear component. Then plot \(E^{(j)}\) versus \(X_j\).

Example

Let’s do a simulation example to see how C+R plots can help distinguish monotone and nonmonotone nonlinearity.

set.seed(1234) # for replication
N <- 100 # set up the sample size
x1 <- rnorm(N, mean = 4, sd = 5)
x2 <- runif(N, min = -5, max = 5)
x3 <- rnorm(N, mean = 10, sd = 2)*(x1/2)
cor(x1, x3)

## [1] 0.9724877

e <- rnorm(N, mean = 0, sd = 1) # set the errors to be normally distributed 
y <- 7 + x1^2 + x2 + 5*x2^2 - 7*x3 + e # set the true y
df <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3, e = e)

In the above chunk, I define the true value \(y\) as a function of \(x_1\), \(x_2\), \(x_3\), and the error term (\(e\)). Specifically, the true regression coefficients are set as follows: \(1\) for the squared transformation of \(x_1\), \(1\) for \(x_2\), and \(-7\) for \(x_3\). Additionally, the term \(5 \cdot x_2^2\) signifies a quadratic effect of \(x_2\), while the intercept is fixed at \(7\). Note, I sepcifically set \(x_3\) to be highly correlated with \(x_1\). To visualize the underlying relationship between \(x_1\), \(x_2\), \(x_3\), and \(y\), let’s start by plotting scatterplots for each pair of \(y\) and \(x\), then running a multivariate regression, and plotting residuals against each \(x\).

fit <- lm(y ~ x1 + x2 + x3, data = df)
df$residual <- fit$residuals

library(tidyverse)

p1 <- ggplot(data = df, aes(x = x1, y = y)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(color = "salmon2", alpha = 0.3) +
  labs(subtitle = "Scatterplot between y and each x") +
  theme_bw()
  
p2 <- ggplot(data = df, aes(x = x1, y = residual)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(method = "lm", color = "steelblue2", alpha = 0.3) +
  labs(subtitle = "Residuals against each x") +
  theme_bw()

p3 <- ggplot(data = df, aes(x = x2, y = y)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(color = "salmon2", alpha = 0.3) +
  theme_bw()
  
p4 <- ggplot(data = df, aes(x = x2, y = residual)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(method = "lm", color = "steelblue2", alpha = 0.3) +
  theme_bw()

p5 <- ggplot(data = df, aes(x = x3, y = y)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(color = "salmon2", alpha = 0.3) +
  theme_bw()
  
p6 <- ggplot(data = df, aes(x = x3, y = residual)) +
  geom_point(color = "gray37", alpha = 0.7) +
  geom_smooth(method = "lm", color = "steelblue2", alpha = 0.3) +
  theme_bw()
  
ggpubr::ggarrange(p1, p2, p3, p4, p5, p6, ncol = 2, nrow = 3)

From the scatterplots of \(y\) against \(x_1\) and \(x_2\), as well as their respective residual plots, it’s evident that while plotting residual plots can indicate whether our model violates the linearity assumption, it doesn’t provide insight into the underlying reasons for the violation. Notably, the residuals plotted against both \(x_1\) and \(x_2\) exhibit a similar convex pattern. Moreover, scatterplots of \(y\) against \(x_1\) and \(x_2\) do not accurately depict their partial relationship. Simply plotting scatterplots between variables fails to consider the influence of other variables in the model, particularly when there are high correlations among predictors. Let’s do the C+R plots instead.

crPlots(fit)

The magenta line represents a smooth curve through the C+R plots, while the blue dashed line represents the regression coefficient of each \(x\) from the multivariate regression model. A notable disparity between the residual line and the component line suggests that the predictor does not exhibit a linear relationship with the response variable. Once again, the advantage of using C+R plots is that they provide information about the partial relationship between \(y\) and each \(x\) while controlling for all other variables. Additionally, they offer insight into how to address nonlinearity. To address the nonlinearity in \(x_1\), applying a square transformation often resolves the issue effectively. As for \(x_2\), we need to consider the quadratic relationship between \(y\) and \(x_2\).

df$x1_square <- (df$x1)^2
fit_update <- lm(y ~ x1_square + poly(x2, 2, raw = TRUE) + x3, data = df)
crPlots(fit_update)

2.2 Violation of Homoskedasticity Assumption

Note: the materials for this section draw from Penn State STAT 501 and here.

The method of ordinary least squares assumes that there is constant variance in the errors (which is called homoscedasticity). The method of weighted least squares can be used when the ordinary least squares assumption of constant variance in the errors is violated (which is called heteroscedasticity). The model under consideration is:

\[Y = X\beta + \epsilon^*\] where \(\epsilon^*\) is assumed to be (multivariate) normally distributed with mean vector 0 and nonconstant variance-covariance matrix.

\[\begin{equation*} \left(\begin{array}{cccc} \sigma^{2}_{1} & 0 & \ldots & 0 \\ 0 & \sigma^{2}_{2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \sigma^{2}_{n} \\ \end{array} \right) \end{equation*}\]

If we define the reciprocal of each variance, \(\sigma_i^2\), as the weights, \(w_i = 1/\sigma_i^2\), then let matrix W be a diagonal matrix containing these weights:

\[\begin{equation*}\textbf{W}=\left( \begin{array}{cccc} w_{1} & 0 & \ldots & 0 \\ 0& w_{2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0& 0 & \ldots & w_{n} \\ \end{array} \right) \end{equation*}\] The weighted least squares estimate is then:

\[\begin{align*} \hat{\beta}_{WLS}&=\arg\min_{\beta}\sum_{i=1}^{n}\epsilon_{i}^{*2}\\ &=(\textbf{X}^{T}\textbf{W}\textbf{X})^{-1}\textbf{X}^{T}\textbf{W}\textbf{Y} \end{align*}\] With this setting, we can make a few observations:

• Since each weight is inversely proportional to the error variance, it reflects the information in that observation. So, an observation with a small error variance has a large weight since it contains relatively more information than an observation with a large error variance (small weight).
• The weights have to be known (or more usually estimated) up to a proportionality constant.
• Weighted least squares estimates of the coefficients will usually be nearly the same as the “ordinary” unweighted estimates. In cases where they differ substantially, the procedure can be iterated until estimated coefficients stabilize (often in no more than one or two iterations); this is called iteratively reweighted least squares.
• The use of weights will (legitimately) impact the widths of statistical intervals.

Example

Let’s simulate the violation of homoskedasticity assumption and apply the method of weighted least squares to address it.

set.seed(1234) # for replication
N <- 1000 # set up the sample size
x1 <- 1:N
x2 <- rnorm(N, mean = 10, sd = 3)
sd <- runif(N, min = 0, max = 2)
e <- rnorm(N, mean = 0, sd = 2*x1) # set the errors to be normally distributed but not at a constant variance
y <- 7 + 3*x1 + x2 + e # set the true y
df <- data.frame(y = y, x1 = x1, x2 = x2, e = e)

In the above chunk, I set the error variance to be more pronounce: as \(x1\) gets larger, the error variance gets larger by multiplying \(x1\). You can see this relationship in the below figure.

ggplot(data = df, aes(x = x1, y = e)) +
  geom_point(color = "gray37", alpha = 0.7) +
  theme_bw()

Let’s now run the regression and test it to see if it violates the homoskedasticity assumption.

fit <- lm(y ~ x1 + x2, data = df) 

	OLS
Variables	Estimates	S.E.	C.I. (95%)
(Intercept)	-65.82	144.57	-349.52 – 217.87
x1	3.13 ***	0.13	2.88 – 3.38
x2	6.73	12.46	-17.71 – 31.17
Observations	1000
R2 / R2 adjusted	0.371 / 0.370
p<0.05   ** p<0.01   *** p<0.001

According to the regression table, you can see that the standard error of \(x2\) is very large, which is due to the heteroskedasticity problem in our data. The estimate of the effect of \(x2\), therefore, is very inefficient.

par(mfrow = c(2, 2))
plot(fit)

Let’s take a look at the third plot. Scale-Location is used to check the homogeneity of variance of the residuals (homoscedasticity). Horizontal line with equally spread points is a good indication of homoscedasticity. This is not the case in our example, where we have a heteroscedasticity problem. It can be seen that the variability (variances) of the residual points increases with the value of the fitted outcome variable, suggesting non-constant variances in the residuals errors.

To formally test for heteroscedasticity, we can perform a Breusch-Pagan test:

# load lmtest package
library(lmtest)

# perform Breusch-Pagan test
bptest(fit)

## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 168.84, df = 2, p-value < 2.2e-16

The Breusch-Pagan test uses the following null and alternative hypotheses:

• Null Hypothesis (\(H_0\)): Homoscedasticity is present (the residuals are distributed with equal variance)
• Alternative Hypothesis (\(H_A\)): Heteroscedasticity is present (the residuals are not distributed with equal variance)

Since the p-value from the test is < 2.2e-16, we will reject the null hypothesis and conclude that heteroscedasticity is a problem in this model.

Since heteroscedasticity is present, we will perform weighted least squares by defining the weights in such a way that the observations with lower variance are given more weight:

# define weights to use
wt <- 1 / lm(abs(fit$residuals) ~ fit$fitted.values)$fitted.values^2

head(lm(abs(fit$residuals) ~ fit$fitted.values)$fitted.values)

##        1        2        3        4        5        6 
## 22.48226 39.03961 48.75519 15.68057 45.26871 47.61134

head(wt)

##            1            2            3            4            5            6 
## 0.0019784279 0.0006561287 0.0004206862 0.0040670185 0.0004879820 0.0004411429

# perform weighted least squares regression
wls_fit <- lm(y ~ x1 + x2, data = df, weights = wt)

	No WLS			WLS
Variables	Estimates	S.E.	C.I. (95%)	Estimates	S.E.	C.I. (95%)
(Intercept)	-65.82	144.57	-349.52 – 217.87	2.14	20.99	-39.04 – 43.32
x1	3.13 ***	0.13	2.88 – 3.38	3.09 ***	0.07	2.94 – 3.23
x2	6.73	12.46	-17.71 – 31.17	0.79	2.55	-4.21 – 5.79
Observations	1000			1000
R2 / R2 adjusted	0.371 / 0.370			0.648 / 0.648
p<0.05   ** p<0.01   *** p<0.001

The weighted least squares model yields a coefficient for \(x_2\) much closer to its true value compared to the original model. Additionally, the residual standard error in the WLS model is smaller than that in the original multivariate linear regression model. These results suggest that the predicted values generated by the weighted least squares model are much closer to the actual observations than those produced by the original model.

The weighted least squares model also has an R-squared of 0.648 compared to 0.371 in the original multivariate linear regression model. This indicates that the weighted least squares model is able to explain more of the variance in y compared to the original multivariate linear regression model.

Discussion IV

TBD

Discussion V

TBD

Discussion VI

TBD